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, ABSTRACT 

Q\ ■ An extension of Schwarzschild's (1979) galaxy-building technique is presented that, 

t-H ' for the first time, enables one to build Schwarzschild models with known distribution 

functions (DFs). The new extension makes it possible to combine a DF that depends 
■ only on classical integrals with orbits that respect non-classical integrals. With such a 

| combination Schwarzschild's orbits are used only to represent the difference between 

. the true galaxy DF and an approximating classical DF. 

' The new method is used to construct a dynamical model of the inner Galaxy. The 

model is based on an orbit library that contains 22 168 regular orbits. The model aims 
7— I ■ to reproduce the three-dimensional mass density of Binney, Gerhard & Spergel (1997), 

£> ' which was obtained through deprojection of the COBE surface photometry, and to 

, reproduce the observed kinematics in three windows - namely Baade's Window with 

(£, b) = (1°, -4°) and two off-axis fields at (8°, 7°) and (12°, 3°). The viewing angle is 
assumed to be 20° to the long axis of the bar and the pattern speed is taken to be 
GOkms^ 1 kpc -1 . 

The model fits essentially all the available data within the innermost 3 kpc. The 
| axis ratio and the morphology of the projected density contours of the COBE bar 

are recovered to excellent accuracy within corotation. The kinematic quantities - the 
' line-of-sight streaming velocity and velocity dispersion, as well as the proper motions 

when available - are recovered, not merely for the fitted fields at (1°, —4°) and (8°, 7°), 
Q ■ but also for three new fields at (8.4°, -6°), (1.21°, -1.67°), and (-1.14°, 1.81°). The 

| dynamical model deviates most from the input density close to the Galactic plane 

c/2 . just outside corotation, where the deprojection of the surface photometry is suspect. 

The dynamical model does not reproduce the kinematics at the most distant window, 
| (12°, 3°), where disk contamination of the data may be severe. 

• ■ . Maps of microlensing optical depth are presented both for randomly chosen stars 

p\ ' and for stars that belong to individual components within the model. While the optical 

?H \ depth to a randomly chosen star in Baade's Window is half what measurements imply, 

the optical depth to stars in a particular component can be as high as the measured 
values. The contributions to the optical depth towards randomly chosen stars from 
lenses in different components are also given. 

Key words: Galaxy: structure - stars: kinematics - Galaxy: Centre - methods: 
numerical 



o 
in 
o 

On 



1 INTRODUCTION 

We live rather far out in the disc of a spiral galaxy, so studies 
of the solar neighbourhood do not provide a balanced view 
of the Galaxy as a whole. It is essential to complement these 
studies with investigations of the inner Galaxy. This more re- 
mote region is hard to study because it is almost completely 
obscured by dust at optical wavelengths. Consequently, the 
first indication that ours is a barred galaxy came from 21- 
cm observations (de Vaucouleurs 1964) . It is only in the last 
few years that a combination of radio-frequency (Binney et 
al. 1991) and near infrared studies (Blitz & Spergel 1991) 



have demonstrated beyond reasonable doubt that the inner 
Galaxy is dominated by a bar whose nearer end lies at pos- 
itive longitudes. Evidence from microlensing has also been 
adduced in favour of a barred inner Galaxy (e.g., Paczyhski 
et al. 1994, Evans 1994), although this remains a controver- 
sial matter (e.g., Bissantz et al. 1997, Sevenster et al. 1998). 

A promising way to constrain the characteristics of the 
bar is to combine observations at various wavelengths with 
theory by trying to construct a dynamical model of the inner 
Galaxy that is compatible with the available data. This is 
the only true test of our assumptions regarding the Galaxy. 
Of course, the current observational data are too scanty to 
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allow us to infer a unique model, and, unsurprisingly, there 
is controversy as to the viewing angle, length and pattern 
speed of the Galactic bar (e.g., Binney et al. 1991, Blitz 

6 Spergel 1991, Sevenster et al. 1998). The real value of a 
dynamical model is that it gives us predictive power and sug- 
gests further observational tests to confirm or constrain the 
structural parameters of the bar. Pioneering studies along 
these lines include those of Pfenniger & Friedli (1991), Sell- 
wood (1993), Zhao (1996) and Sevenster et al. (1998). 

In this paper we build a dynamical model by populating 
orbits in a given potential. This technique was pioneered by 
Schwarzschild (1979, 1982, 1993) and further developed by 
Merritt & Fridman (1996) and Zhao (1996). Schwarzschild's 
original application of his method was to the building of 
models of elliptical galaxies in which most of phase space was 
regular. In general, rotating barred potentials support very 
considerable numbers of irregular orbits and a straightfor- 
ward application of Schwarzschild's method is not profitable. 
In Section 2 we therefore extend Schwarzschild's method in 
such a way that both regular and chaotic regions of phase 
space can be populated. 

Our model of the inner Galaxy is constrained to fit 
the three dimensional luminosity density of Binney, Ger- 
hard & Spergel (1997), which was designed to reproduce 
the infrared data from the COBE satellite. It was recovered 
by Richardson-Lucy deconvolution of the COBE data, after 
correction for extinction using a model in which dust is dis- 
tributed throughout the Galaxy (Spergel, Malhotra & Blitz 
1996). The resulting three dimensional density is specified 
on 59 x 59 x 39 data-cube, corresponding to a box that is 
lOkpc on a side and 2.8 kpc thick. Our dynamical model 
is also constrained to reproduce the kinematic observations 
in three windows - Baade's Window with (I, b) — (1°, —4°) 
and the two off-axis fields of Minniti et al. (1992) at (8°, 7°) 
and (12°, 3°). All our constraints are described in Section 3, 
where particular attention is paid to the best ways of com- 
paring models to kinematic data. To have value, a dynam- 
ical model requires a large number of orbits, which are the 
basic building blocks in Schwarzschild's method. Our orbit 
library, which contains the density and kinematic contribu- 
tions of over 23 000 orbits, is described in Section 4. Section 
5 discusses the penalty and merit functions used to drive the 
mass on the orbits towards the desired input density. Our 
final dynamical model of the inner Galaxy is analyzed in 
detail in Section 6. There we predict the values of measur- 
able kinematic quantities in several fields near the Galactic 
centre, describe the model's DF, and analyze the variation 
of optical depth to gravitational microlensing both with po- 
sition on the sky and with the Galactic component to which 
either the source star or the lensing object belongs. Section 

7 sums up. 



2 EXTENDING SCHWARZSCHILD'S 
TECHNIQUE 

2.1 Splitting the DF 

In Schwarzschild's original technique, we calculate N orbits 
labelled by i = 1, . . . , N in the given potential and determine 
the fraction pij of the time that the ith orbit lies in each 
of the j = 1, . . . , K cells. Let p° bs be the system's original 



density in the jth cell with volume Vj. Then we seek the 
non-negative weights u>; that minimize the discrepancies 



A, = V 3 p° hB - 



M 
~N 



N 

y^ W i Pij- 



(1) 



Here, the constant factors can be absorbed into the defini- 
tion of the weights, but we have written them explicitly for 
later convenience. 

So long as the galaxy model is stationary in inertial 
space or in a frame of reference rotating at constant pat- 
tern speed, the Hamiltonian H is an isolating integral of the 
equations of stellar motion. If the system is axisymmetric, 
one component of angular momentum, L z , will also be an 
isolating integral. Consequently, by Jeans' theorem, any DF 
that is a function of H, and where appropriate of L z , will 
satisfy the collisionless Boltzmann equation (see e.g., Bin- 
ney & Tremaine 1987). This suggests that we approximate 
the true DF by a classical DF, / class (#) or / class (#, L z ) 
as appropriate. Since H and L z are known functions of the 
conventional phase-space coordinates 



(r,p), 



(2) 



we can readily calculate the value of the density or any kine- 
matic quantity to which J class would give rise. Except in 
special cases, the true DF, /, will depend on non-classical 
isolating integrals. Therefore we write 



f(w) = r iass (H(w),L z (w)) + r*(w), 



(3) 



where / rcg is a function that depends on non-classical iso- 
lating integrals. We show below how Schwarzschild's tech- 
nique can be used to determine simultaneously J class and 
/ rcg under the assumption that non-classical isolating inte- 
grals exist only along orbits that are regular in the sense 
that they have three effective isolating integrals. The latter 
we take to be actions Ji, i = 1,2,3. For simplicity of expo- 
sition we henceforth assume that the only classical integral 
is H, which is certainly the case for the inner Galaxy. Since 
DFs that depend on H alone generate isotropic velocity dis- 
tributions in the appropriate frame of reference, we refer to 
the component generated by J class as the isotropic compo- 
nent and call its DF the isotropic DF, / iso = / class (#). All 
the results below generalize trivially to more general classi- 
cal DFs. 

Note that in our convention f(w) is a probability den- 
sity and is always normalized to unit phase-space integral 
inside the box used to fit the density (see below). 



2.2 From weights to DF 

The key problem involved in representing / rcg = / — / lso 
with Schwarzschild's technique, is the determination of the 
value of the DF that is implied by a given set of weights w™ 15 . 
We solve this problem by using an arbitrary, everywhere 
positive, correctly normalized probability density f s (w) to 
sample points in phase-space. With these points as initial 
conditions, we integrate along orbits. Then we define 



f(j(w)) 



lim -|- / 

Jo 



lim - / dtf (w(t)) 



if the orbit 
is regular, 

otherwise, 



(4) 
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where w(t) is the orbit with initial conditions w. Con- 
structed thus, f s (J) automatically obeys the collisionless 
Boltzmann equation. Hence, any distribution function con- 
sisting of regular orbits only can be written as 



r s (j) 



(•») 



where w ICS (J) is some weight-function. Moreover, f s (J) 
gives the relative probability that a regular orbit obtained 
by sampling phase space according to the probability den- 
sity f B (w) will be the orbit with actions J . The correspond- 
ing probability-density for picking the ith regular orbit, is 
(2nf{N B /N Teg lp(J i ), where iV s initial conditions gave rise 
to jV reg regular orbits, while the factor of (2tv) 3 accounts for 
the phase-space volume at constant J. 

The isotropic part of the DF, / ls °, is most conveniently 
represented by a superposition of basis functions Bi(Ej), for 
example splines, 



(6) 



such that the complete DF is determined by the set of 
weights K°, ra rcE (J)}. 



2.3 The non-negativity constraint 

In irregular regions of phase space, / rcg vanishes and the 
total DF, /, is equal to the isotropic DF, / lso . Hence the 
latter must be non-negative in irregular regions. By equa- 
tion (|^), the requirement that the total DF be non- negative 
translates to the following constraint on the weights: 

f iso (H(J)) y . w) so B, (H(J)) 



f s (J) 



f s (J) 



In general, / ISO will be everywhere non-negative and this 
equation allows w rcg (J) a limited degree of negativity. 

The physical implication of an orbit having a negative 
weight is this. Frequently, both regular and irregular orbits 
exist at a given value of H . In a three-dimensional system, ir- 
regular orbits tend to Arnold diffuse over much of the phase- 
space hypersurface of constant H. They are, however, rig- 
orously excluded from regions of the hypersurface that are 
occupied by regular orbits. The isotropic part of the DF may 
assign non-zero density to irregular orbits at some value of 
H when the total phase-space density, /, in some regular 
region of the same energy is negligible. In such a case the 
weights of the regular orbits would approximately satisfy the 
equality condition in equation ((?]) . Hence, regular orbits can 
be important even if no stars are on them by virtue of their 
ability to exclude stars on irregular orbits from subsets of 
their energy hypersurfaces. 



2.4 From observations to weights 



To calculate /, the weights {wf 



1 (J)} must be deter- 



mined from observational data. Any observable moment, 
such as the density, is linear in the DF, and can be split 
into the contributions arising from / lso and / rcg . That is 



where H(w) is the function that characterizes the moment 
in question - see below. The left-hand side of equation (ph 
is given by observations and the right-hand side is a linear 
function of the weights. A sufficiently large set of such equa- 
tions, for different functions n, can now be treated as an 
inverse problem and a standard technique used to recover 
the weights {iu| so , w rcg ( J)} from the measured values of the 
left-hand sides. 

Since H(w) is a known function, the mom ents n[/ lso ] 
are straightforward to evaluate (see Section 4.2). Using 
and the standard identity d 6 w = d 3 Jd 3 #, where the 6i are 
the angle variables conjugate to the Ji, the moments arising 
from / rog can be written as 



(9) 



n[f cg ] = J d 3 j«/ cg (j) J d 3 en(w(j,e)). 

The right part of the right-hand side can be calculated nu- 
merically via the time-averaging theorem (e.g., Binney & 
Tremaine 1987): 

TT(j) = ( 27T )~ 3 J d 30 n(«0 = T lim -f J dt 6(w(<)). (10) 

The integral over d 3 J in equation (^) is performed by Monte- 
Carlo integration employing N Tes regular orbits sampled 
from the probability density to (2ir) 3 (N s /N Icg ) f s (J) that 
was computed above. This yields: 



(ii) 



where uf* s is the weight of the ith orbit and is the time 
average of H(w) on this orbit. Note that in practice we do 
not know w ICS (J) but wj eg , and hence ui rcg (un(i)) , for the 
sampled orbits only. 

As an example of this formalism, consider M J rog , the 
contribution from regular orbits to the mass in cell j. In 
this case we define 

1 if to lies in cell j, 
otherwise. 



(12) 



Then the phase-space integral over n^, i.e. the left-hand side 
of equation (^), becomes the fraction of the total mass that 
is in cell j. Here, "total mass" refers only to the mass Ma 
inside the box that contains all cells - recall that f(w) is nor- 
malized to unit phase-space integral inside that box. Hence, 
multiplying equation ([ll]) by Ma we have the mass Mj cg on 
regular orbits in cell j 

JV r °8 

Mrcg Ma \ rcg rcg 

i=0 

where 

P :; g = lim i / diHj ( Wi (t)) 
T ^°° 1 Jo 

is the fraction of the time that the ith regular orbit is in the 
jth cell. The analogy between this equation and equation 
(nl) is clear. 



(13) 



(14) 



n[/] 



d 6 ™ /(«,) n(w) = n[/ is °] + n[f 



(8) 
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Figure 1. This shows the geometry of the bar in the Galactic 
plane. The three dimensional luminosity density of Binney et al. 
(1997) resides in a box that is 5 kpc square in the plane. The 
viewing angle of the Sun is taken as 20° . The x-axis corresponds 
to the bar's major axis, the y-axis is the minor. The intercepts 
of the line of sight through the Galactic Centre with the box are 
marked. The bar is rotating in the clockwise direction. 



3 APPLICATION TO THE GALAXY 

We will be working in a rotating frame of reference, in which 
the dynamics are governed by the Hamiltonian 



H(w) = |p 2 + $(r) - O, ■ (r X p), 



(15) 



where p is the canonical momentum per unit mass and ft is 
the angular velocity of the bar. Standard manipulations en- 
able us to recast this into the form (e.g., Binney & Tremaine 
1987) 

HO) = |« 2 + $ cff (r), $ cff (r) = $(r)-i0 2 (x 2 + y 2 )(16) 

where v is the velocity in the rotating frame and $ c fi is the 
effective potential. The Hamiltonian H , which is an exactly 
conserved quantity, is also known as Jacobi's integral Ej. 



3.1 The density 

As the original density p°, we use the three-dimensional 
model that Binney et al. (1997) obtained by deproject- 
ing the near-infrared COBE surface photometry. The sur- 
face photometry employed by Binney et al. had been cor- 
rected for absorption by dust by a procedure that is outlined 
in Spergel, Malhotra & Blitz (1996). The non-parametric 
Richardson-Lucy algorithm used by Binney et al. is based 
on the assumption that the density is eight-fold symmetric 
(with respect to the major, intermediate and minor axes). 
While this assumption is reasonable enough for the bar it- 



self, it causes features like spiral arms to be incorrectly re- 
produced - see Binney et al. for a discussion of this prob- 
lem. Both the viewing angle and the pattern speed of the 
Galactic bar are somewhat controversial (e.g., Binney et al. 
1991, Sevenster et al. 1998). For the sake of definiteness, we 
take the bar's viewing angle as 20° and the pattern speed 
as 60kms -1 kpc -1 . The Sun is assumed to lie at a Galac- 
tocentric radius of 8 kpc and 0.014 kpc above the Galactic 
plane. The circular speed at Ro is taken to be 200 kms -1 . 
The Sun's peculiar motion is (10, 5, 7) kms -1 , where the first 
component is along the line of sight towards the Galactic 
centre, the second is in the direction of Galactic rotation 
and the third component points towards the north Galactic 
pole. In other words, the Sun is moving in towards the Galac- 
tic centre, leads the local standard of rest and moves up and 
away from the Galactic plane (see e.g., Binney & Merrifield 
1998). All this means that in the rest frame of the Galaxy 
with (x, y, z) coordinates aligned with the symmetry axes of 
the bar, the Sun has phase-space coordinates 



(x,y,z) = (7.5,2.7,0.014) kpc, 
{ Px ,p y ,Pz) = (60.7,-196,6.98) kms" 1 . 



(17) 



The other coordinate system we frequently use is a helio- 
centric frame. The line-of-sight distance s and Galactic lon- 
gitude and latitude (I, b) are the configuration-space coor- 
dinates. Velocity space is given by radial or line-of-sight ve- 
locity v\ os , together with the proper motions (fj,i, p,j,). In this 
system, the Galactic centre has coordinates 

(s,£, b) = (8, 0, 0) kpc, 

wioB = -10kms" 1 I (18) 
(Hi,fJ,b) = (-5.4, -0.18) masyr -1 

Note that at the Galactic centre, 1 masyr -1 =37.92 kms -1 . 

The density plays a dual role in Schwarzschild's method: 
it both constrains the weights through equation ( ^c| ) and, 
through Poisson's equation, specifies the potential in which 
the orbits are calculated. The Galactic density of Binney 
et al. (1997) comes in two parts. First, inside a box that 
is 10 kpc on a side and 2.8 kpc thick and whose geometry 
and position w.r.t. the Sun is shown in Fig. ^, the density 
is specified on 59 x 59 x 39 grid points. This density distri- 
bution was obtained by non-parametric Richardson-Lucy 
deprojection of the photometry, followed by multiplication 
by a constant mass to light ratio T. Between grid points, the 
density is evaluated through three-dimensional cubic splines. 
Outside the box the density is given by an analytic func- 
tion that is based on the work of Spergel, Malhotra & Blitz 
(1996): 

P(r) = Po[fb(r) + U(r)] (19) 
where 

2 / 2 



h = fo 



(l + a/ao) 1 - 8 ' 



Zll 



Zl 
2 \ 1/2 



7?de -fl/fld , 



a = ( x 2 + ^ + — ) and R= ^ x 2 + 



(20a) 
(20b) 

(20c) 
these 



The constants in 

equations are as follows: po = 2.05 x 1O 8 M0 kpc -3 for the 
L-band, /o=624, a m = 1.9kpc, ao = 100pc, i?d = 2.5kpc, 
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Figure 2. The left panel shows contours of the effective potential in the Galactic plane. The right panel shows the effective potential 
along the major, intermediate and minor axes. The minimum of the effective potential is at the Galactic Centre, while maxima occur 
both on the j/-axis and at (x, y) = (±2.9, ±2.2) kpc. There are saddle points on the a;-axis and close to the maxima along the y-axis. 



z = 210 pc, 2i = 42 pc, a = 0.27, f] = 0.5, and £ = 0.6. Physi- 
cally, these equations specify an exponential disk, in which 
the vertical distribution is the sum of two exponentials, and 
a triaxial bulge that extends to R ~ 2 kpc. Since we use the 
density distribution only at R > 5 kpc, the bulge has a neg- 
ligible impact on our model, box and the density exterior to 
the box is achieved by setting p = ppbox + (1 — <?)pcxt, where 



exp 



< s < 1 



(21) 



and s 4 = (x/5kpc) 4 + (j//5kpc) 4 + (z/1.4kpc) . Appendix 
A describes the technique used to determine the gravita- 
tional potential <3> generated by this mass distribution. The 
interesting features of the potential are conveniently de- 
scribed in the rotating (r,v = r) frame of the bar. Contours 
of the effective potential (j^i]) are plotted in Fig. ^. The plot 
is dominated by a basin whose rim lies at i?.~3.8kpc. The 
lowest point on the rim defines the critical value of the po- 
tential, 3?™- Orbits whose Jacobi energy Ej is smaller •J?™ 
cannot cross the rim. Consequently, orbits belonging to the 
isotropic component with Ej < ^eff 1 are °f two types: those 
that lie entirely inside the rim, and those that lie outside it. 
The latter extend to infinity, as do orbits with Ej > $3f\ 
Since such unbound orbits are useless for galaxy modelling, 
/ ISO is non-zero only for Ej < Q^fp, and then describes orbits 
that lie entirely within the rim. 

In Fig. ^| a maximum centred on (2.9, 2.2) kpc is con- 
spicuous. The potential's other stationary points are the 
Lagrange points Li to L5 that were first identified in the 
context of the restricted three body problem. Li_2,3 lie along 
the i-axis, with Li being at the Galactic Centre and L/2,3 
opposite each other at x~ ± 3.7 kpc, while L.4,5 lie on the 
y-axis at y w ± 3.5 kpc. 



3.2 The kinematics 

3.2.1 The selection function 

The main problem when observing stars in the bulge is ob- 
scuration. This problem is worst within the Galactic plane, 
and a relatively unobscured optical view at the Galactic cen- 
tre is only possible in a very few windows, the most famous 
of which is Baade's Window at (£, b) = (1°, -4°). 

To compare observational data with a dynamical model, 
careful thought has to be given to the selection criteria 
applied to obtain the sample of observed stars. The selec- 
tion function e(s, M) is the probability that a star of abso- 
lute magnitude M that lies at heliocentric distance s is in- 
cluded in the sample. Unhappily, for many published surveys 
e(s, M) is hard to determine. If a survey contains all stars 
brighter than the limiting magnitude m max and fainter than 
some cutoff magnitude some magnitude m m i n , the function 
e(M, s) is given by 



i{M,s) = 



< M + 5 log yjj^ + -ys <m 

max 

: otherwise 



(22) 



where 7 is the differential extinction in magnitudes per unit 
distance. Direct observations of stars and their colour excess 
yield the total extinction A in certain windows, but 7 itself 
is only available from three-dimensional dust-models whose 
reliability is uncertain. Even if a dust model fits the dust dis- 
tribution well when the latter is averaged over some scale, 
it will be inaccurate at individual points because extinction 
within the disk is very patchy. Therefore, in equation ( pj[ ) we 
use the measured total extinction A rather than the differ- 
ential extinction 7 of a model - since most extinction lies in 
a foreground screen close to the Sun, this procedure should 



© 0000 RAS, MNRAS 000, 000-000 



6 Ralf M. Hdfner, N. Wyn Evans, Walter Dehnen and James Binney 



not give rise to large errors. To allow for the patchiness of 
extinction, we treat A as a Gaussian random variable. 

If the <f>(M) is the stellar luminosity function, then the 
general selection function e(s) 



10 o 



dM <p(M) e(M,s), 



(23) 



gives the probability that a star at distance s that has un- 
known absolute magnitude will included in a survey. Com- 
bining the last two equations, we obtain the family of generic 
selection functions that we will use: 



dAe 



-A 2 /2a 



-A-51og[s/10pc] 



(24) 



dM <j>(M). 



-A-51og[s/10pc] 



The samples in the bulge are normally dominated by giant 
stars. Each type of star (e.g., K giants or M giants) may 
be considered to have a relatively narrow band of intrinsic 
luminosities M that is well-modelled by a Gaussian 



4>{M) 



exp 



(M - M int ) 2 



(25) 



Here, the average intrinsic luminosity Mint and dispersion 
om vary according to stellar type. 

3.2.2 The fitted windows 

The kinematic data are fitted in three windows. These are 
Baade's Window and the two off-axis windows studied by 
Minniti et al. (1992) at {£, b) = (8°, 7°) and (12°, 3°). These 
windows are chosen because the selection criteria seem rea- 
sonably clear-cut and reproducible. 

Baade's Window has been studied extensively in the 
visible wave-band and it is the only one for which proper- 
motion data are available. Sharpies, Walker & Cropper 
(1990) measured line-of-sight velocities for an unbiased sam- 
ple of 239 late-type M giants. The stars were divided into 
two groups which show different kinematics. The first group 
contains 14 bright stars (/ < 11.8) with a relatively small 
velocity dispersion of 71 jli'j' kms -1 . They are believed to be 
either in the outer part of the bulge on the solar side, or 
foreground disk giants, or younger, more massive asymp- 
totic giant-branch stars. The second group, of fainter stars, 
is attributed to the bulge itself. Sharpies et al. argue that 
this group can be considered complete. The velocity dis- 
persion is considerably higher at HSl^kms -1 . The mean 
velocity is 4±8kms _1 . The value of 113lgkms _1 for the 
line-of-sight velocity dispersion seems to be very robust as 
other studies find similar values (Rich 1988). This dataset 
is modelled with one of the generic selection functions in- 
troduced in equation (|24|) with parameters J max = 13.4 and 
/mm = 11-8. The mean extinction Ai = 0.87 and its disper- 
sion o"a =0.1 are estimated from Table 3.21 in Binney & 
Merrifield (1988). 

Spaenhauer, Jones & Whitford (1992) conducted the 
first, and so far only, survey of stars in Baade's Window 
that measured proper motions rather than just line-of-sight 
velocities. They compared plates taken at epochs that are 20 
years apart and selected 800 stars with B — V > 1.4. These 




Figure 3. Top left: A set of selection functions defined by 
m m j, n = m max - 1, Mnt = 0, A = 0, a M = 0.5 and <t a = 0.1. 
?"max is taken to be 13, 14, 15 and 16, respectively. Top right, 
bottom left, bottom right: Two panels are shown for each 
of the windows at (1° -4°), (8°, 7°), and (12°, 3°), respectively. 
The upper panel shows the probability of finding a star at dis- 
tance s in a given sample while the panel on lower panel shows 
the density along the line of sight. The sharp cut-off occurs when 
the edge of the box is reached. 



Table 1. The parameters used in the section functions for the 
data-sets provided by Sharpies et al. (1990), Spaenhauer et al. 
(1992) and Minniti et al. (1992). 



(i,b) 




^min 


A a A 


"lint 




Ref 


(l°,-4°) 


13.4 


11.8 


0.87 0.1 


-2.75 


1.0 


Sh 


(l°,-4°) 


17.5 


16 


1.8 0.2 


1.3 


1.0 


Sp 


(8°, 7°) 


16.5 


12 


1.5 0.5 


-1.5 


1.0 


M 


(12°, 3°) 


16.5 


12 


1.5 0.5 


-1.5 


1.0 


M 



The first column gives the Galactic coordinates of the window. 
The next two columns report the faint m ma x and bright m m i n cut- 
offs used in the selection function. The fourth and fifth columns 
give the mean extinction A and its dispersion a a - The mean in- 
trinsic magnitude Mi n t and its dispersion <7jYf f° r the luminosity 
function of the sample are reported in columns six and seven. 
These refer to the wave-bands in which the observations are taken. 
The last column is the mnemonic for the window: Sh for Sharpies 
et al., Sp for Spaenhauer et al. and M for Minniti et al. 



are almost entirely K and M giants. Out of these 800 stars, 
371 were excluded due to overcrowding or because they were 
too faint. Stars between B = 17 and B = 19 (roughly) are in- 
cluded. Spaenhauer et al.'s sampling seems best-reproduced 
by making the following assumptions about the generic se- 
lection function: V max = 17.5, V m i n = 16, A v = 1.8, a a = 0.2, 
Mv = 1.3 and om = 1-0. It should be noted that neither 
Vniax nor Vmin are hard limits but rather soft boundaries 
derived from the final sample, while the average extinction 
is taken from Stanek (1996). Spaenhauer et al. used their 
measured stars to define the reference frame and so their 
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means in the proper motions are intrinsic. Only the disper- 
sions in the proper motions - namely ai = 3.2 ± 0.1 masyr -1 
and Ub = 2.8 ±0.1 mas yr~ - carry dynamical meaning. 

Finally, Minniti et al. (1992) studied stars in two off- 
axis fields, namely (I, b) = (8°, 7°) and (12°, 3°). Stars with 
R < 16.5 were pre-selected. Out of these only the reddest 
were taken. The limits in colour were chosen to correspond to 
the locus of K stars. They estimate that in the two fields ~ 10 
and ~ 30 per cent are disk stars. The peculiar motion of the 
Sun is assumed to be 15.4 kms -1 towards (£, b) = (51°, 23°). 



The results are vi oa — 45 ± 10 km s 



, = 85±7kms _1 for 



the first field, and wi os = 77 ± 9 kms~ , <ti os = 68 ± 6 kms 
for the second field. When comparing the predictions of our 
model with the Minniti et al. data, the following generic 
selection function is used: i? max = 16.5, i? m i n = 12, Ar — 1.5 
and a a = 0.5. The mean intrinsic magnitude of the K stars 
is taken as Mr = — 1.5. The mean extinction Ar is merely a 
crude estimate and is open to debate. The parameters used 
in the generic selection functions are listed in Table |l| for all 
the three windows. 

As an illustration of the generic selection functions, we 
plot an example in Fig. |[ These selection functions are of 
one magnitude width. In other words, m max = m m i n — 1. The 
top left panel shows the selection functions plotted against 
heliocentric distance. The top right, bottom left and bottom 
right panels are for the three windows (1°, —4°), (8°, 7°), and 
(12°, 3°), respectively. In each case, the upper figure shows 
the number of stars n(s) picked up at a heliocentric distance 
s, namely 



n(s) oc e(s) p(s) s 



(26) 



and the lower figure shows the density along the line of sight. 
The important point is that the selection function has a 
crucial influence on the observed kinematics, as it controls 
where the stars are picked up. As m m i n is varied from 13 to 
16, the observables are dominated by foreground stars, then 
bulge stars proper, and finally stars lying behind the bulge. 
In principle, the selection function enables us to probe kine- 
matic structure along lines of sight. Instead of just a single 
value of a variable such as the line-of-sight velocity disper- 
sion, there is really a function depending on the magnitude 
cut-off. 



4 THE ORBIT LIBRARY 

This section is concerned with the numerical implementation 
of our extension of the Schwarzschild algorithm. The first 
sub-section discusses the choice of the sampling distribution 
function, which controls the selection of the regular orbits in 
the orbit library. The second sub-section discusses the choice 
of the isotropic component. 

In this and subsequent sections, the energy E or the 
Hamiltonian H(w) — Ej are given in units of GMq kpc -1 ~ 
(978 kms -1 ) 2 , while the unit of the angular momentum is 
VGMq kpc ~ 978 kpc km s -1 with G denoting Newton's 
constant of gravity. 



4.1 The sampling distribution function 

Since we use the isotropic DF to populate the irregular parts 
of phase space, our orbit library contains only regular orbits. 



Figure 4. Histograms of the distribution of normalized Liapunov 
exponents for orbits within three narrow ranges of values of the 
Jacobi energy Ej. The total area under each of the curves is 
normalized to unity. The dotted vertical lines show our division 
into regular and irregular orbits at these Jacobi energies. Orbits 
to the left of the dotted line are regular, orbits to the right are 
irregular. 



Regular orbits in a rotating barred potential fill a relatively 
small volume of phase space around the closed prograde or- 
bits and a rather larger volume around the closed retrograde 
orbits. We ensure that the sampling distribution f s is large 
in these regions by using a sum of products of a density 
distribution p(r) of an ellipsoidal Hernquist model, together 
with functions h(p) of momentum only: 



f=Y,A iP {r) hi{p). 



(27) 



Here, the non-negative numbers Ai may be chosen for con- 
venience subject to the condition 1 = J d 6 w f. 

In the first component of the sampling distribution / s , 
hi(p) is strongly peaked around the momentum of closed 
prograde orbits. For the second component, /i2(p) is peaked 
around the velocity of closed retrograde orbits. For the third 
component, hz{p) has a broad peak around p = 0. We refer 
to these components as the prograde, retrograde and hot 
components, respectively. 

Mathematically, for the hot component, h is given by 



(2tt) 3 / 2 



■ exp 



Px 

2<j1 



P 2 y 



2al 



Pz 

2al 



(28) 



where the ai are specified by Table ^. 

In our coordinate system, the Galactic bar and disk 
have negative angular momentum (see Fig. ^). We define 
p,j, to be p, (which is not the momentum conjugate to 
azimuth <j> but the Cartesian momentum p resolved in the 
direction of increasing <f>) and assume that on prograde or- 
bits p ~ —v c e^, while on closed retrograde orbits p = v c e r j > , 
where v c is defined by 



v c (R) = 0.25 



1 + 



^ O.lkpc y 



(29) 



For the prograde and retrograde components, h is given by 



hi,: 



1 



(2^)3/2^(7^(7, 



exp 



A. 

2*4 



2al 



Pz 

2ol 



(30) 



where the plus sign is taken for the prograde component, 
and the minus sign for the retrograde component. Table ti 
ives the values of the parameters that appear in equation 



The spatial parts of the sampling density (^) are de- 
fined by 



p(r) 



with 



10 1.5 kpc 

27-7T m 



i 1 + 1.5 kpc) 



m 2 = x 2 +y 2 + (z/0.4) 2 . 



(31) 



(32) 



Orbits are followed for ~ 200 dynamical times. Over this 
time, the Jacobi energy is conserved to typically one part in 
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Table 2. The velocity dispersions in kms 1 of the three compo- 
nents of the sampling distribution function / s (p7[). 

Component crj? cta <t x cr y a z 



which does not exist in a barred galaxy, but is denned to be 
(in the units given above) 



prograde 

retrograde 

hot 



60 
60 



60 
60 



100 100 



50 
50 
60 



The dispersions, given in kms -1 , are used in the sampling func- 
tion and provide the basis for choosing the orbits in the prograde, 
retrograde and hot components of the orbit library. 

10 6 . Liapunov exponents A are used to distinguish regular 
from irregular orbits on the principle that A = for a regu- 
lar orbit. The process of estimating A is discussed by Udry 
& Pfenniger (1988). In practice, we extrapolate to infinite 
time by fitting the estimate X(t) obtained by following the 
orbit for time t to A = A + b/t. These unnormalized Liapunov 
exponents are converted to normalized ones by multiplying 
by the orbital time. This is defined to be the mean time 
between successive passages through the plane y = 0. Fig. ^ 
shows three histograms of values of normalized Liapunov ex- 
ponents Anorm for orbits in three narrow ranges of Ej. In all 
three cases, there is a sharp peak around Anorm = 0, which 
corresponds to the regular orbits. After integrating for an in- 
finite time, one might expect all irregular orbits at any one 
Jacobi energy Ej to be equivalent. Since we only calculate 
for finite time, the peak corresponding to the irregular orbits 
is broadened. The dotted vertical lines in Fig. ^ separate the 
regular and the irregular orbits. In addition to this, we also 
include all orbits whose Liapunov time is greater than five 
bar rotation times. This seems reasonable as such orbits do 
not evolve on the time-scale (~ 100 bar rotation periods) of 
interest to us. 

The prograde sampling function is used to select 200 000 
initial conditions. These give rise to only 5 713 regular orbits. 
Similarly, 100 000 initial conditions selected from the sam- 
pling distribution of the hot component give rise to 5 512 
regular orbits. Only 50 000 initial conditions selected from 
the sampling function of the retrograde component give rise 
to 10 943 regular orbits. The final orbit library contains 
jV rcg = 22 168 regular orbits. The library records the prob- 
abilities of each orbit being in any given cell and the time- 
averaged sampling density. 

Fig. |^ shows the density of the library's orbits in a con- 
venient projection of orbit space. In a conventional Lindblad 
diagram for an axisymmetric galaxy, the angular momen- 
tum L z of orbits is plotted horizontally, and their energy 
is plotted vertically. The top panel of Fig. [B] is a modified 
Lindblad plot for our system, in which orbit-averaged val- 
ues of L z and energy are plotted horizontally and vertically. 
These averages are effective rather than classical integrals. 
Since the Hamiltonian satisfies H — E — flL z , the classical 
integral H is a linear combination, (E) — Q{L X ) of the effec- 
tive integrals. We choose to plot (E) rather than H because 
the former is a truer guide to an orbit's physical size than 
the latter. In a classical Lindblad diagram, the allowed re- 
gion \L Z \ <L C (E), where L C (E) is the angular momentum 
of a circular orbit of energy E, tapers as one descends to 
smaller values of E. In Fig. H this tapering has been largely 
suppressed by plotting horizontally not (L z ) but {L z )/L c , 
where L c is not the angular momentum of a circular orbit, 



M<£>) = 6oo(^ + - 25 ) 



(33) 



Since the Galactic bar has a negative pattern speed, 
prograde orbits lie on the left hand side of Fig. [| and con- 
tours of constant H slope from top left to bottom right. The 
red curve in Fig. |B| is the contour for H — $™. The region 
of strong chaos associated with corotation is evident in the 
white sea that cuts the (E) axis between —0.035 and —0.075. 
Within this sea there is a long thin island of regularity. The 
panels marked CI and C2 show two orbits within this island. 
One is trapped around the standard maximum of <& e fT at L4. 
The other is trapped around the non-standard maximum of 
$ e ff at (±2.9, ±2.2) kpc. At smaller values of (E), a ridge of 
orbits trapped around the prograde, bar-supporting x\ fam- 
ily is apparent near the left-hand edge of Fig. B . The panels 
marked A, Bl and B2 in the lower half of Fig. | show repre- 
sentative orbits from this region. Panel D shows a prograde 
orbit at the largest value of (E) plotted in Fig. ^. This is a 
typical disk orbit. 

On the right-hand, retrograde side of orbit space, orbits 
of the X4 family occupy a region of regularity that extends, 
unbroken, from the smallest to the largest energies. The pan- 
els labelled K, L and M show representative orbits from this 
region. 

In Fig. ^| colour shows the density of orbits, while white 
dots show the orbits themselves. Many sharp chains of dots 
can be discerned: these mark the paths of stable resonances. 



4.2 The isotropic component 

The isotropic DF can be conveniently represented as a linear 
combination of some basis functions - the coefficients in this 
expansion are the free parameters that then determine / lso . 
We have employed second-order basis-splines (see e.g., Stoer 
& Bulirsch 1980). Hence, 



AT lso 

rw = ^»rB,(HW), 

=1 



(34) 



where N' B ° = 1000 and 



Bi(Ej) = 



AE~i 



Ej — -E.Ti-i for Ej g \E H _ U E H \, 
Ej l+1 - Ej for Ej g [Ej u Ej i+1 ], (35) 
otherwise. 
Here, AEj denotes the grid spacing, taken to be constant, 
while the are constants that enforce some chosen normal- 
ization. The basis functions Bi(H(w)^J can be interpreted 
as building blocks containing all orbits with Jacobi integral 
Ej = H(w) for which Bi(Ej) > 0. 

As mentioned in Section 3.1, outside corotation, irregu- 
lar orbits can escape to infinity. Since these irregular orbits 
are occupied via / ls °, we must restrict / ls ° to values of Ej 
of orbits that cannot cross the potential rim at corotation. 
This means that (i) / lso is non-zero only for Ej < ^"g 1 , and 
(ii) at fixed Ej < "I?™ refers only to orbits inside corotation. 
Thus, strictly speaking equation (^) for / lso contains two 
Heaviside functions, one ensuring Ej < 5™ and the other 
:r|<x r i m , where 3>™ = ^ e s(xtim, 0, 0). For simplicity, we 
have suppressed these Heaviside functions. 
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Figure 5. Top panel: The density of orbits in the library of regular orbits is colour-coded, with black the lowest density. Scaled, 
time-averaged angular momentum is plotted horizontally and time-averaged energy vertically. The red curve is the contour H = 3>S? , 
above which orbits are not confined by the Hamiltonian and the isotropic part of the DF, / lso , is set to zero. Individual orbits are marked 
by white dots. Lower panels: representative orbits. The location of each orbit in the top panel is indicated by the letter to the right or 
the orbit and the ((E), (L z )) coordinates on top of it. Orbits are followed for 200 dynamical times. 



Table 3. Labelling of the velocity moments in the merit function. 



in cell j at any time. The weighting of the cells, pnormj, 
chosen to be 



1 . . .3 <«i os ) , (fit) , (fj, b ) 
4... 6 (*, los 2 ) , (a e 2 ) , ( M6 2 ) 

7 ... 9 (v loB U e ) I (flosMb) , (PlPb) 

Not all these components of the velocity dispersion tensor are 
available for all the windows. The most complete information is 
known for Baade's Window — namely, v\ os , <ti os , erg and u\y. 



Computing moments of the isotropic DF, / lso , reduces 
to computing them for the Bi. For the mass density, for 
example, 



p ls °(r) = M n ^ wl i°Pi ( r )i 

i 

where 



p?°(r)= d 3 vBi(Ej 



15 AEj 



(36) 



[DU-2D!-D 5 i+1 ] (37) 



with £>| = En — 4> c ff (r)- The ki are determined by normal- 
izing pf° such that 



1=1 d 3 r pT{r)- 



(38) 



This integral over real space is done numerically by the 
Monte-Carlo method. 



5 THE MERIT AND PENALTY FUNCTIONS 

The combined mass density of the orbits must match the 
original density p° bs used to create the potential (see Section 
3.1). To this end we minimize the merit function 



obs 

Pi 



(41) 



where N SBm j is the number of building blocks contributing 
to grid cell j. The effect of this weighting scheme is a min- 
imization of relative errors in the density modified by the 
Poisson error expected from the sampling. 

The kinematics of our model are given by the matrix 
iu'jl with i = 0, ... ,9 representing the ith velocity moment in 
window j, as defined in Table |§| Appendix B gives details 
of the calculation of the moments. The contribution to the 
merit function of the kinematic constraints is 



=EE^ 



(42) 



where 7^ is a weighting matrix and the superscripts m and 
o refer to the model and observed quantities, respectively. 
In some windows only a subset of the moments is available. 
For all three of the constraint windows (described in Section 
3.2), line-of-sight velocity vi os and dispersion ai os are known. 
Only in Baade's Window are the dispersions at and a;, also 
available. 

The problem posed by the minimization of the merit 
function is very ill-conditioned, so regularization is neces- 
sary. This is achieved by enforcing a smoothness constraint 
- 'neighbouring' orbits should have 'similar' weights. A suit- 
able penalty function is the mean-square value of the sec- 
ond derivative of the logarithm of the weights with respect 
to some distance - this vanishes when the dependence of 
weights on distance is a power law. For orbits belonging to 
the isotropic DF, a suitable distance is provided by Ej, and 
the penalty function reads 



E 



where the sum extends over all cells and 



M° 



i=l 



ISO ISO 



JV r °s 

-y- 



,rcg rcg 



(39) 



(40) 



where Mp=5.18x 1O 1O M , while M° bs is the mass con- 
tained in the jth cell of our mass model according to BGS. 
pf° is the integral of pf°(r) [equation over cell j. Be- 

cause the Bi are normalized to unit phase-space integral, 
pf° is just the probability for a star whose phase-space co- 
ordinates are drawn from Bi(H(w)^J to be found in cell j 
at any time. Equivalently, pjj° s , defined by equation (Q, is 
the probability that a star on the ith regular orbit is found 



TV" 



E 



In wf 



21n«4 so + In wj+i 



(A£j)= 



(43) 



For the regular orbits, the distance measure is provided by a 
set of effective integrals (c.f. Merritt & Fridman 1996, Zhao 
1996). Specifically, we employ the time averages I\ = (p 2 z ), 
l2 = {L z ) and 13 = (E) - note that 73 is the mean orbital 
energy, not the value of Ej. The penalty function for regular 
orbits then reads 



preg 



V 1- 

E 



ln W r g -<ln» rcg ) A(l 



(d) 



Mi) 



(44) 



Here, A(i) is a neighbourhood of orbit i and denotes 
the average within this neighbourhood. A(i) is defined to 



© 0000 RAS, MNRAS 000, 000-000 



10 Rolf M. Hdfner, N. Wyn Evans, Walter Dehnen and James Binney 



u 
ex 



> 2 r 




o 

Q_ 



N 



2 3 
x [kpc] 



4 




1 


2 




1 










8 







6 







4 







2 





2 3 4 
x [kpc] 




1 2 3 

y [kpc] 



4 



Figure 6. Logarithmic contours of the projected mass in the principal planes of the Binney et al. (1997) model (dotted lines) are 
compared with those of our dynamical model (full lines). The units are such that the total mass in the 10 kpc X 10 kpc X 2.8 kpc box is 
unity. 



consist of the 8 orbits closest to orbit i. The distance, dij, 
between orbits i ^ j is defined by 

3 r v, , 

(45) 



d 



E 

k=l 



where er/ fc is the dispersion of Ik over all the regular orbits. 

The final quantity to be minimized, Q, is a linear com- 
bination of the density and kinematic merit functions, Q den 
and Q kin , and the penalty functions P Icg and P iso . The fi- 
nal numerical factors in this linear combination were chosen 
as follows. The relative weight for Q km was chosen to be 
as small as possible without increasing the deviation from 
the observed kinematics by more than the observational un- 
certainty (1<t). Similarly, the relative weight for the penalty 
functions were chosen as large as possible without worsening 
the fit to the mass density by more than 2 per cent overall. 

To enforce the non-negativity of /, we substitute the 
weights uij SO and u>J eg by fij and /3j in the following way: 



^ 1st 
: In Wj 



In 



(46a) 



(46b) 



Hence, /3 r ' 



tends to the lowest value com- 



patible with the positivity constraint i 

qreg-i 



Minimizing Q with 

respect to {/3j S °, /3^° s } always results in a physical model. 
Unfortunately, two attractive properties of the original op- 
timization problem are lost in the transition from the w's 
to the f3's: Q dcn no longer depends linearly on the variables, 
and the boundedness of the solution is not guaranteed (as 
the (3's can diverge). 

Due to the great number of orbit weights to be deter- 
mined, the memory requirement of the adopted optimiza- 



tion algorithm must not increase with the number of un- 
knowns faster than linearly. This excludes fast schemes such 
as non- negative least square fitting (e.g., Zhao 1996). An- 
other popular choice, the iterative Richardson-Lucy method 
(Newton & Binney, 1984; Statler, 1987), is not applicable to 
our problem because its kernel is not positive definite. Our 
final choice fell on the conjugate gradient algorithm (e.g., 
Stoer & Bulirsch 1980; Press et al. 1988). It satisfies the 
stringent memory requirements and is an improvement on 
the steepest descent method in so far as the directions in 
which it does its line minimizations are conjugate to each 
other. 



6 RESULTS 

This section describes our dynamical model of the Milky 
Way's bar. Section 6.1 describes how the model reproduces 
the density and kinematic constraints. Section 6.2 provides 
kinematic predictions of the model, while Section 6.3 dis- 
cusses the phase-space structure. 



6.1 The constraints 

6.1.1 The density 

The full contours in Fig. ^ show the density of our final 
dynamical model projected onto the three principal planes. 
The broken contours show the corresponding projections of 
the input density of Binney et al. (1997). Overall, the fit be- 
tween the two density distributions is good, with the average 
discrepancy ~ 10%. The biggest contribution to this error 
comes from the innermost two layers of cells close to the 
Galactic plane, and particularly around corotation on the 
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Figure 8. The mass profiles of our model are shown along the x, y and the z-axes. The isotropic DF contributes two-thirds as much 
mass as the regular component near the centre. Along the rr-axis the contribution from the isotropic DF decreases rapidly, while along 
the z-axis the isotropic component becomes steadily more important until it dies out near z ~ 1.2. Relatively few regular orbits pass 
through the minor axis, so the discrepancies between our model and the input density are substantial further than 1.2 kpc down the 
z-axis. 



CJ 



> 2 




Figure 7. The mass within 0.14 kpc of the Galactic plane is 
shown for the Binney et al. (1997) input density (dotted con- 
tours) and for the present model (full contours). While the model 
can reproduce the mass in the inner parts very well it fails to re- 
produce the over-densities on the y-axis that are most likely due 
to symmetrized spiral arms. 



y-axis. Outside these areas, the error is very much smaller. 
The model reproduces both the axis ratio and the detailed 
shapes of the contours in the (x,z) and (y,z) planes very 
well out to the limits of the box, although there are some 
small discrepancies on the z-axis - see below. 

Fig. [?] highlights the discrepancies between the model 
and input densities in the (x, y) plane by contouring for each 
model the mass within 0.14 kpc of the plane. The two sets 
of contours agree well within corotation, but significant dis- 



crepancies occur further out. These discrepancies are not 
worrying for two reasons. First, obscuration is at its worst in 
the Galactic plane, and in the region one cannot confidently 
deproject the surface photometry. Second, the deprojection 
algorithm of Binney et al. assumes that all structures are 
eight-fold symmetric with respect to the principal planes of 
the bar. Spiral arms are less symmetric and spurious features 
will arise when one attempts to deproject a spiral distribu- 
tion under the assumption of eight-fold symmetry. Indeed, 
Binney et al. show that features remarkably like the 'ob- 
served' the density maxima along the j/-axis arise when the 
Binney et al. algorithm is used to deprojection a four-armed 
spiral distribution. Englmaier & Gerhard (1998) provide fur- 
ther evidence for this interpretation by showing that better 
fits to the observed longitude-velocity diagrams for HI and 
CO are obtained when the flow of gas in the inner Galaxy 
is calculated in a potential that includes contributions from 
the density maxima. Moreover, the contributions from the 
maxima induce a four-armed rather than a predominantly 
two-armed spiral in the gas. Hence several lines of evidence 
indicate that the density maxima in the deprojected disk are 
artifacts associated with four spiral arms. Such spiral arms 
inevitably lie beyond the reach of our modelling technique. 

Fig. |^ shows density profiles of the model along the 
x- and z-axes. At the centre the isotropic component con- 
tributes over 30 per cent of the mass. Within the plane this 
component diminishes rapidly in importance as one moves 
away from the centre. Along the z-axis, by contrast, the 
isotropic component gains in importance as one moves from 
the centre because many of the orbits passing through the 
z-axis are irregular. Beyond ~ 1.2 kpc, the isotropic compo- 
nent is unable to contribute to the density because orbits 
inevitably have Jacobi energies in excess of "3?™) and the 
discrepancies between our model and the deprojected pho- 
tometry increase. 



© 0000 RAS, MNRAS 000, 000-000 



12 Rolf M. Hafner, N. Wyn Evans, Walter Dehnen and James Binney 

Table 5. Model predictions for unmeasured quantities. 
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(1°, -4°) 




0.26 


-0.01 




3.5 


2.8 


-0.05 


0.03 


0.04 


Sh 


(1° -4°) 


3.8 


0.25 


-0.01 








-0.05 


0.03 


0.04 


Sp 


(8° 7°) 




-0.85 


-0.22 




6.0 


3.4 


0.52 


0.05 


-0.14 


M 


(12° 3°) 




0.70 


-0.10 




4.4 


2.1 


0.08 


-0.05 


0.01 


M 


(8.4° -6°) 


77 


0.87 


0.29 


78 


3.1 


3.0 


0.31 


-0.01 


0.15 


TT 


(1.21°, -1.67°) 


8 


-0.16 


-0.09 


134 


3.5 


3.0 


-0.11 


0.02 


0.04 


Bl 


(-1.14° 1.81°) 


-8 


-0.20 


0.01 


112 


4.4 


2.6 


-0.10 


-0.01 


0.03 


Bl 



Units are masyr -1 for proper- motions and kms - 1 for velocities. The dj are the dimensionlcss 
correlation coefficients of the observable velocity dispersion tensor, that are defined by equation ([47]). 
The first four windows are those used to constrain the model (data used for that purpose are given 
in Table Q and are omitted here), the fourth window has been observed by Tiede & Tendrup (1997), 
and the last two by Blum et al. (1997). The line of sight streaming for the Blum et al. fields has not 
been corrected for the reflex motion of the Sun to enable comparison with his data. 



Table 4. Kinematic data used to constrain the model. 





Quantity 


observed 


Model 


Ref 


(1°, -4°) 




4±8 


4 


Sh 






113±5 


115 




(1° -4°) 


°"los 


120 


111 


Sp 




°l 


3.2±0.1 


3.6 








2.8±0.1 


2.8 




(8° 7°) 


"los 


45±10 


45 


M 




°los 


85±7 


80 




(12°, 3°) 


%>s 


77±7 


75 


M 




°"los 


68±6 


95 





The units are kms" 1 for velocity dispersions and masyr" 1 for 
proper-motion dispersion (cr^ and 0-5). The last column is the 
mnemonic for the observers (see Table |l|). The biggest discrepan- 
cies occur at (12°, 3°), where our dynamical model has too much 
dispersion and not enough streaming. 

6.1.2 The kinematics 

It is essential to apply both photometric and kinematic con- 
straints to obtain plausible models of the inner galaxy. The 
dynamical model is required to reproduce the data provided 
by Sharpies et al. (1990) and Spaenhauer et al. (1992) in 
Baade's Window, together with the observations of Minniti 
et al. (1992) at (8°, 7°) and (12°, 3°). Note that we main- 
tain a distinction between the Sharpies et al. (1990) and the 
Spaenhauer et al. (1992) data-sets, as the model is viewed 
at the same window through different selection functions. 
Table ^] shows how the model fares. In Baade's Window, the 
line-of-sight streaming velocity, «i os , and velocity dispersion, 
<ti os , together with one of the proper- motion dispersions <ib 
are well reproduced. The remaining proper-motion disper- 
sion a i is higher than measured by Spaenhauer et al. (1992). 
This is because our orbit library, and hence our model, prob- 
ably contains too many retrograde orbits. Turning to the 
Minniti et al. (1992) fields, the streaming velocity and the 
dispersions at (8°, 7°) are reproduced to within the error 
bars. The more distant window at (12°, 3°) is more difficult 
to get right. One worry is that disk contamination is likely 
to be severe in this outer window, which may mean that 
Minniti et al.'s (1992) results need correction. This explana- 
tion is consistent with the fact that our model has a higher 



dispersion and a lower streaming velocity than are suggested 
by the data. Overall, though, Table ^ encourages us in the 
belief that our dynamical model is a good representation 
of the inner Milky Way and that it is useful to make some 
kinematic predictions from the model. 

6.2 Kinematic predictions 

Let us stay for the moment with our constraint fields. In 
Baade's Window, the means in the proper motions pn and 
fib are not given by Spaenhauer et al. (1992), but it may 
become possible to recover them at some time in the fu- 
ture. The mixed components of the velocity dispersion ten- 
sor are also thus far unmeasured, although there is prelimi- 
nary claim of a measurement of the vertex deviation from a 
small sample by Zhao, Spergel & Rich (1994). In the Minniti 
et al. (1992) fields, only line-of-sight quantities are available. 
Table @ presents the predictions of our model for all the un- 
measured quantities. [We have not given results for (12°, 3°) 
because our dynamical model does not reproduce the exist- 
ing data there.] We also compute the correlation coefficients 
between the observable velocity dispersions 

222 

„ _ a los,< „ _ °"los,6 r , _ cr eb . . 

O s f = , Csb = , Cib = • (4/) 

Zhao, Spergel & Rich (1994) suggest that C s e is a useful di- 
agnostic of bulge triaxiality. At Baade's Window, this quan- 
tity vanishes for a steady-state axisymmetric density distri- 
bution. The predictions for these correlations are presented 
in the tables. They are, of course, measures of the misalign- 
ment of the principal axes of the velocity dispersion tensor 
with the (s,£, b) axis set. In Baade's Window, the dispersion 
tensor has principal semi-axes in the ratio 0.83 : 1 : 0.77, 
with the longest axis pointing almost in the ee direction. 
Hence, Baade's Window is a rather poor place to look for 
the signature of triaxiality. Our barred model has a disper- 
sion tensor whose alignment is almost the same as that of 
an oblate axisymmetric model! In the Minniti et al. (8°, 7°) 
field, the dispersion tensor is very strongly anisotropic, with 
semi-axes in the ratio 0.25 : 1 : 0.50. The long axis points al- 
most in the ei direction and has the high value of 238 kms" 1 . 
In the model, there are more and more retrograde stars 
picked up as one moves further from the Galactic Centre. 
The velocity dispersion in the longitudinal direction rises as 
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the prograde and retrograde stars become present in almost 
equal numbers. 

Now, let us see how our model fares in comparison 
with data for two new fields. Tiede & Terndrup (1997) 
present the results of a study of 189 stars in a field at 
{£, b) = (8.4°, —6.0°). The selection criteria are described in 
detail in their paper but are not simple to reproduce. A crude 
approximation to their selection procedure is to use the 
generic selection function (^) with parameters m max = 17.0 
and m m i n = 12. Tiede & Terndrup provide values for the 
mean extinction A =1.1 and its dispersion aA =0.2. They 
measure the line-of-sight dispersion cri os of their sample 
to be 75±lkms~ 1 . Table ^ shows the predictions of our 
model. Good reason for believing in the model's reliabil- 
ity within corotation is that this value of the dispersion is 
very reproduced. The diagonalised tensor has semi-axis ra- 
tios 0.42 : 1 : 092 and is rather strongly misaligned with the 
(s,£,b) axis set, as the large correlation C s t indicates. 

Blum et al. (1994) studied stars in two fields very 
close to the Galactic Centre at (1.21°, —1.67°) and at 
(—1.14°, 1.81°). The stars comprising the sample were se- 
lected in a manner that tried to exclude disk stars. It is not 
so easy to reproduce their selection procedure. A crude ap- 
proximation is to take the selection function as unity for all 
heliocentric distances s satisfying 6 kpc < s < 10 kpc. From 
their data, Blum et al. deduce that wi OB = 14 ± 23 km s _1 and 
cti os = 128± 14kms _1 at (1.21°, -1.67°) and that v los = - 
75±24kms" 1 and a los = 127 ± 17 at (-1.14°, 1.81°). The 
predictions of our model for these two new fields are reported 
in Table [|. Again, there is the reassuring circumstance that 
both the dispersions are reproduced to within the errors. 
One of the streaming velocities is recovered to within the 
error bars, but one is not. At (1.21°, —1.67°) the disper- 
sion tensor has semi-axes (141, 126, 112) kms -1 , while at 
(-1.14°, 1.81°) the semi-axes are (111, 167, 100) kms" 1 . In 
both cases the principal axes are not strongly misaligned 
with the (s, I, b) coordinate directions, so, as in Baade's Win- 
dow, triaxiality will be hard to establish unambiguously in 
these fields. 

Finally, Figs. | to [t(] show the kind of kinematic data 
that may become available in the very near future. Here, 
we have imagined that our dynamical model is observed 
through four windows using the generic selection function 
( pi| ) with a width of one magnitude. We assume that the 
extinction and intrinsic magnitude of the stellar population 
vanish, so m max is in effect the maximum distance modulus 
in the sample. It is apparent from Figs. ^|to |l^that a wealth 
of additional information is uncovered when the kinematics 
is studied as a function of apparent magnitude. 

As one application of the figures, let us examine how to 
improve the slight discrepancies between the model and the 
observations in the constraint fields. The quantity at is too 
large by a factor of ~ 5 per cent in Baade's Window. From 
Fig. ^, it is evident that this can be corrected by slightly 
increasing the faint cut-off, as the curve of at versus m max 
falls for m max > 14.5. If, for example, either the total extinc- 
tion or the faint cut-off is less severe than we have assumed, 
then at is lowered to give better agreement with the Spaen- 
hauer et al. data. The quantity at is also slightly too large in 
our dynamical model, but Fig. ^ makes it clear that increas- 
ing the faint magnitude cut-off makes little difference to its 
value. Some of the properties of the curves in Figs. ^| and [H] 



(8.4", 6") 




Figure 9. The left panel gives the kinematic predictions at 
Baade's Window (1°, —4°) and the right panel at the field stud- 
ied by Tiede & Terndrup (8.4°, —6°). Here, the generic selection 
function is used with m m ; n =m E 



1. 



(8',7-) 



(12\3°) 




Figure 10. Kinematic predictions at Minniti et al.'s windows. 
The left panel shows results for (8°, 7°), the right panel for 
(12°, 3°). Again, the generic selection function is used with 



are readily explained. The line-of-sight streaming v\ OB typi- 
cally gives a one-humped curve, as it is greatest when stars 
are picked up by the selection function at roughly the tan- 
gent point. The line-of-sight dispersion a\ OB can sometimes 
give complex curves, as in the (12°, 3°) field of Minniti et 
al. Here, as we move along the line of sight, the selection 
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Figure 11. The distribution function of the regular component of 
the final model is shown projected into the (L z )/L c — (E) plane. 
The phase-space density is largest in the bottom left corner of the 
figure, which is occupied by small, nearly harmonic box orbits. At 
higher energies the density peaks along the ridge of the prograde 
ici-orbit family. 



Figure 12. The distribution of the regular component's mass 
within the (L z )/L c — (E) plane. The colour of each cell encodes 
the sum of the weights assigned to the orbits that lie in that cell. 



function picks up first a mixture of prograde and retrograde 
stars in the near side of the disk, then mainly prograde stars 
in the bar proper, and finally prograde and retrograde stars 
in the far side of the disk. This causes a\ os to rise, then fall, 
and then rise again. 

Clearly, as sample sizes increase, it will become increas- 
ingly inappropriate to characterize the kinematics of a given 
field by just mean velocities and dispersions. Even observa- 
tions of external galaxies can now deliver other measures of 
the line-of-sight velocity distribution (LOSVD), such as the 
Gauss-Hermite coefficients /13 and h^. Still richer informa- 
tion should be available in the case of the Milky Way be- 
cause, by varying m m i n and m max , we can probe the LOSVD 
at different points along the line of sight. More sophisticated 
data sets will surely resolve much of the degeneracy that now 
plagues dynamical models. 



6.3 The distribution function 

Fig. [n] shows the distribution function of the regular com- 
ponent of the final model projected onto a plane whose axes 
are averaged energy (E) and averaged angular momentum 
(L z ), normalized by the reference value, L c , that is defined 
by equation ( |33| ) . The phase-space density is highest at bot- 
tom left, and decreases with increasing energy for both pro- 
grade and retrograde orbits. These trends are in the same 
sense as the variation of / ISO , which decreases monotonically 
as H increases from bottom left to top right. 

Fig. [12] shows the distribution of the mass of the reg- 
ular component within the ({L z }, (E)) plane: the colour of 
each pixel encodes the sum of the weights u) rcg assigned to 
orbits whose effective integrals place them within that pixel. 
[Equation @ shows that the mass on an orbit is propor- 
tional to its weight.] Comparing Figs, [n] and [li] we see that 
mass is concentrated at much larger energies than phase- 
space density. This is a simple reflection of the fact that the 
amount of phase-space volume that is associated with a pixel 
in Figs. [0] and [13 increases rapidly with {E} . The surpris- 
ing feature of Fig. [l2| is large concentration of mass in the 
top right corner of the diagram. This mass is on retrograde 
orbits around corotation. It lies there because at these ef- 
fective energies the isotropic component cannot contribute, 
and there is a lack of regular prograde orbits. 

We find that 75 per cent of the mass in the box is on 
orbits that are confined to be inside corotation. Of these, 
about a sixth (12 per cent) are in the isotropic component, 
half (37 per cent) are prograde regular and a third (26 per 
cent of the total in the box) are retrograde regular. Thus 
Within corotation, it is prograde orbits that are dominant. 



Table 6. Mass and optical depth towards Baade's Window (BW) 
and the Galactic Centre (GC). 



component 



Sources by comp. 
GC BW 



Lenses by comp. 
GC BW 



isotropic 
prograde 
retrograde 
hot 

total 



0.12 
0.32 
0.39 
0.17 

1.00 



1.1 10~ 5 1.1 10~ 6 

2.3 10~ 5 1.3 10~ 6 

4.6 10~ 5 2.0 10~ 6 

1.3 10~ 5 1.3 10~ 6 



3.5 10~ 6 3.6 10- 7 

2.9 10" 6 5.1 10~ 7 

3.910" 6 2.710" 7 

9.410~ 6 1.1 10~ 7 

2.0 10~ 5 1.210 



-6 



The prograde, retrograde, and hot components refer to the reg- 
ular orbits with (L Z )/((L Z — (-£<z)) 2 } 1/ ' 2 smaller than —3, larger 
than 3, or in between. The relative masses only refer to the con- 
tributions inside our box of 10 X 10 X 2.8 kpc. The left columns 
give for the sources drawn from the various components, while 
the lenses are taken from the entire model (inside the box). The 
rightmost two columns are for sources throughout the model and 
lenses in individual components. 



We have tried and failed to build models that place a much 
larger fractions of mass on the isotropic component. Inside 
corotation, chaotic orbits play quite a small role, although, 
as Fig. |i| illustrates, they are indispensable along the mi- 
nor axis. Outside corotation, retrograde orbits dominate the 
mass budget. The second column of Table ^ gives the mass 
contained in each component. 

Fig. ^ gives an overview of the distribution of mass be- 
tween the various components by plotting AM/AE, the mass 
per unit increment in energy. Since the DF of the isotropic 
component is a function f lso (H) of the Hamiltonian, cal- 
culating the corresponding form of AM/AE involves deter- 
mining the volume in phase space within which both E(w) 
and H(w) lie within specified ranges, and then multiplying 
this volume by / lso and integrating over all H. Details of 
this calculation are given in Appendix The long-dashed 
curve in Fig. [l3| shows the resulting curve, which rises fairly 
steadily with E until it drops sharply to zero as the cutoff 
at H = $"ff n cuts in. 

The full curve in Fig. ^| shows AM/AE for the regular 
component in the approximation that each orbit contributes 
mass only to the energy that is equal to (E), the time- 
average of E(w) along the orbit. Since E generally does not 
vary greatly along a regular orbit, this is a good approxima- 
tion. Comparing the full and long-dashed curves in Fig. |l3| 
we see that inside corotation the regular and isotropic com- 
ponents make comparable contributions to the overall mass, 
with the isotropic component dominant at the lowest ener- 
gies, and the regular component mostly dominant further 
out. The short-dashed and dotted curves in Fig. [l^ show 
the contributions to the regular component from prograde 
and retrograde orbits, respectively. Prograde orbits generally 
contribute more than half the mass. Around the corotation 
energy, E — —0.05, there is a glaring exception to this rule, 
however, as the contributions to AM from prograde orbits 
plunges to zero and retrograde orbits provide all the mass. 
At slightly higher energies the relative importance of the two 
orbit types reverses, sharply. 
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Figure 13. The differential mass distribution of the isotropic and the regular component. The regular part is further split into prograde 
and retrograde orbits. The isotropic component stops rather abruptly at E = —0.09 since it is confined to inside corotation. 



full mode 




Figure 14. Logarithmic (to the base 10) contours of the microlensing optical depth for /3 = 0. Left: lenses drawn from the full model 
while the source population differs from panel to panel; this corresponds to the results on the left part of Table ^j. Right: sources drawn 
from the full model while the lens population differs from panel to panel. Note that the density, and hence the optical depth, of the 
component with the isotropic DF ends abruptly at corotation. 



6.4 The microlensing optical depth 

The optical depth for gravitational microlensing, r, is given 
by 



[l[?a D °P*(?*) f dDd D4D s ~D d )p d (D d ) 



AtvG o 



D. 



2/3 



(48) 



dD, 



D 2 B p s (D a 



D, 



id 



Here, the subscripts d and s refer to deflector and source 
objects, respectively and /3 parameterizes the efficiency with 
which stars are picked up along the line of sight (see e.g., 
Kiraga & Paczyriski 1994). The value /3 = is appropriate 
for sources luminous enough to be included in a microlensing 
survey no matter how far away they lie within the Galaxy, 
while 13 = 1 is appropriate for less luminous sources that have 
a probability of being included in the survey that falls off 
with heliocentric distance s as s~ . Typically, (3 — Q is ap- 
propriate for red-clump stars, while (5 ~ 1 is appropriate for 
main-sequence stars. The top left panel of Fig. [l4| shows 
contours of optical depth for the case (3 — 0, while the 
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bottom entries in the last two columns of Table ^ give nu- 
merical values for the direction towards the Galactic centre 
and Baade's Window. These depths are only for microlens- 
ing a source that lies within the box by a lens that also lies 
within the box. Hence, before they are compared with obser- 
vations, these numbers should be augmented by the contri- 
bution from the foreground disk, which for Baade's Window 
typically amounts to 0.1 - 0.3 x 10~ 6 (e.g., Evans 1995). The 
augmented value for Baade's Window, namely 1.4 x 10~~ 6 , is 
smaller than the observational value: Udalski et al. (1994) 
found t — (3.3 ±1.2) x 10 -6 for clump stars, while Alcock et 
al. (1997) found r = (3.9±i;|) x 10" 6 for clump stars near to 
Baade's Window. This problem was discussed in detail by 
Bissantz et al. (1997), who examined essentially the input 
density distribution. 

The four right-hand panels of Fig. 14| and the other 
entries in the last two columns of Table qbreak the over- 
all optical depth down into contributions from individual 
components. Towards the Galactic centre, the largest con- 
tribution comes from lenses in the hot component, while in 
Baade's Window this component makes the smallest contri- 
bution because of its high degree of central concentration. 
Nearly half the overall optical depth in Baade's Window 
comes from lenses in the prograde component. 

The four lower left panels in Fig. [H] and the second and 
third columns of Table H give the optical depth values when 
the lenses are drawn from the full model, but the sources are 
separated by component. The table shows that these vary 
by a factor of up to 2 in a given direction. Specifically, in 
Baade's window, sources that move on retrograde orbits are 
nearly twice as likely to be lensed as sources that belong 
to the isotropic component, and, in fact, have an optical 
depth that is compatible with the measurement of Udalski 
et al. The origin of this phenomenon is simple: objects that 
lie near the outside of the box, but diametrically opposite 
to the Sun, inevitably have larger optical depths than ob- 
jects at the Galactic centre. Fig. |l^ shows that the largest 
radii are dominated by retrograde orbits, while the isotropic 
component is concentrated at the centre. There is a real pos- 
sibility that the density of clump stars does not rise towards 
the Galactic centre as rapidly as the overall stellar mass den- 
sity, with the result that their optical depth is anomalously 
high like that of the model's retrograde objects. 

Barred models can possess asymmetric microlensing 
maps (Evans 1994, Zhao & Mao 1996). Two competing ef- 
fects can be discerned in Fig. [l4|. First, lines of sight to 
sources at negative longitudes are longer and pass through 
more of the dense central bulge than lines of sight at posi- 
tive longitudes. This is because the sources at negative lon- 
gitudes are on average further away. Second, the selection 
effect, controlled by the parameter /3, can enhance the op- 
tical depth at positive longitudes rather than negative, as 
here the sources appear to be brighter. These two effects 
tend to work against each other, so the final asymmetry of 
the map is a subtle matter. For example, in the panel on 
the extreme right of the middle row of Fig. |l4| which is for 
lenses that lie in the prograde component and sources that 
are drawn from the full dynamical model, we see enhanced 
optical depths at negative longitudes, but only outside the 
plane. Within the plane, any interpretation is difficult be- 
cause this is where our dynamical model is least reliable. By 
contrast, when the sources lie in the retrograde component 



near the Galactic plane (extreme left panel of bottom row of 
Fig. ^) , the asymmetry goes the other way. Here, the main 
contribution is from sources at the far end of the Galaxy 
outside corotation, and now the optical depth is enhanced if 
the lenses lie roughly halfway between observer and source. 
This favours positive longitudes, as there is a ready supply 
of such lenses in the foreground bar. 

From the top-left panel of Fig. jl4] we see that the optical 
depth is always larger at negative longitudes when |6| > 1°. 
The coming decade is likely to see the gradient of microlens- 
ing optical depth with Galactic longitude and latitude mea- 
sured, at least at certain spots in the Bulge. For example, 
the third panel from the left in the bottom row of Fig. |l4| 
shows that when the lenses are drawn from the retrograde 
component, the gradient with respect to Galactic longitude 
near the plane is negligible, although the latitudinal gradient 
remains steep. At on-axis (6~ 0°) fields with £ > 10°, the op- 
tical depth is mainly provided by the retrograde component. 
Hence, gradient information, by giving an indication of the 
shape of the contours in the microlensing maps, may pro- 
vide clues as to which components are providing the largest 
numbers of sources and lenses. 



7 CONCLUSIONS 

We have shown how Schwarzschild's method can be used to 
produce models with known DFs. Previously, models con- 
structed by Schwarzschild's (1979) method have had known 
DFs only when the DF depended only on classical integrals 
(e.g., Cretton et al. 1999). The great merit of Schwarzschild's 
method is precisely its ability to handle the general case in 
which non-classical integrals are important, so our algorithm 
for determining the DF of a general model must be counted 
a significant advance. 

One advantage of being able to determine the DF of a 
Schwarzschild model is that one can then combine a crude 
approximation to the galaxy with a DF that depends only 
on classical integrals, with a more general DF obtained 
by Schwarzschild's method. In this approach, only the dif- 
ference between the true DF and the classical DF need 
be reproduced with orbits. The resulting model will have 
higher resolution in real- and velocity-space, and be easier 
to interpret physically than one constructed by the classical 
Schwarzschild technique. 

A key point is that when a classical DF is used in the 
construction of a Schwarzschild model, the weights of orbits 
can be negative because the orbits may be subtracted from 
the underlying classical component, subject only to the con- 
straint that the total phase-space density is non-negative. 
When regular orbits are assigned negative weights, these 
orbits are important because they are excluding stars on 
chaotic orbits from certain regular parts of phase space. 

We have illustrated these new techniques by using them 
to construct a model of the central kiloparsecs of the Milky 
Way. This problem is less well suited to the new techniques 
than the classical problems of modelling axisymmetric sys- 
tems and triaxial systems with negligible figure rotation 
because there is only one classical integral, the Jacobi en- 
ergy, and this integral is not confining for the more ener- 
getic stars. Moreover, the phase space of the inner Galaxy 
is more chaotic than regular, and Schwarzschild's technique 
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works best in a highly regular phase space. From the fact 
that we have succeeded in constructing a reasonable model 
of the inner Galaxy notwithstanding these difficulties, we 
infer that the new techniques would solve easier classical 
problems with some facility. 

Our new Galaxy model reproduces essentially all of 
the available density and kinematic data within corotation 
(~3.6kpc). In particular, the three-dimensional bar density 
of Binney et al. (1997) is faithfully reproduced in the inner 
parts. At corotation, our dynamical model does not repro- 
duce the over-densities on the minor axis found by Binney 
et al. - but these are probably artifacts of the deprojec- 
tion algorithm, which assumes that the underlying model is 
eightfold symmetric and must misrepresent spiral features. 
The kinematic data in Baade's Window (1°, —4°) and in the 
field studied by Minniti et al. (1992) at (8°, 7°) are all fitted 
to within the observational uncertainties. The model does 
not fit data for the outer window of Minniti et al. (1992) at 
(12°, 3°) accurately in that it has a higher line-of-sight dis- 
persion and a lower streaming velocity than the data. This 
mismatch may be caused by disk contamination in the Min- 
niti et al. sample or by too low an assumed value for the 
mean extinction in our model. 

In addition to fitting the initially prescribed data, the 
model furnishes many predictions. We have tested a few of 
these predictions against data in the literature and find that 
they are consistent with measurements by Blum et al. (1994) 
at (1.21°, -1.67°) and (-1.14°, 1.81°), and by Tiede & Tern- 
drup (1997) at (8.4°, -6°). For several well studied fields we 
predict values for quantities, such as proper-motion disper- 
sions, that have not yet been measured. We find that in 
central fields, such as Baade's Window, the principal axes 
of the velocity-dispersion tensor tend to be approximately 
aligned with the line of sight. Hence, these are not good 
locations to probe for kinematic signatures of triaxiality. 

The most puzzling aspect of our model is its reliance on 
retrograde orbits in a narrow band of energies, that corre- 
sponds to radii around 3kpc. This fact has an observation- 
ally testable consequence: at (8°, 7°) the dispersion tensor 
should be extremely anisotropic, with its longest axis aligned 
nearly with the ee direction, because samples should contain 
roughly equal numbers of prograde and retrograde stars, on 
roughly circular orbits. Is this counter-intuitive prediction 
an artifact of our model or a robust prediction? It is a con- 
sequence of two inputs: (i) the density profile of Binney et 
al., and (ii) the decision to exclude from the orbit library 
orbits with Liapunov times shorter than five bar rotation 
periods. 

The density distribution of Binney et al. is very uncer- 
tain at the relevant radii by virtue of a combination of the 
effects of spiral structure and obscuration. Photometry of 
external barred galaxies, such as NGC 1300 (Elmegreen et 
al. 1996) suggests that the Binney et al. profile may be sig- 
nificantly in error, and that the density should drop steeply 
inside corotation to a lower, approximately level value that 
extends throughout the highly chaotic region around coro- 
tation. An approximately constant density around corota- 
tion could be provided by prograde orbits that have apocen- 
tres well outside corotation. By specifying a density gradient 
around corotation, we have obliged the fitting algorithm to 
employ orbits with apocentres at smaller radii, and the only 
available regular orbits are retrograde. Alternatively, we may 



have been too conservative in our criterion for excluding or- 
bits as chaotic: an orbit that has a short Liapunov time may 
never the less remain trapped for a substantial fraction of a 
bar rotation period. In a subsequent study it would be inter- 
esting to use the spectral approach to determining orbital 
regularity (Binney & Spergel 1984; Carpintero & Aguilar, 
1998) instead of Liapunov exponents, which are expensive 
to compute and may be more misleading in cases of mild 
stochasticity. Moreover, even though half of all disk galaxies 
have bars, it is not clear that any individual bar lives for 
a Hubble time; bars may dissolve and reform on a shorter 
timescale. 

We have calculated microlensing maps for our model. 
These show that the overall microlensing optical depth of 
luminous sources in Baade's window is at least a factor 2 
smaller than the observations require. They show also that 
sources in different components have optical depths that dif- 
fer by up to a factor 2. Hence, it is important to characterize 
the populations to which stars that are known to have been 
lensed belong, and to understand the mix of populations 
that characterizes the stars that are regularly monitored for 
microlensing events. In principle, our model allows one to de- 
termine the distribution of durations of the lensing events in 
any field, and we plan to present such distributions shortly. 

It would be useful to extend the work in this paper 
by building dynamical models for different bar morpholo- 
gies, viewing angles and pattern speeds. However, our ex- 
periments have convinced us that there is very considerable 
freedom to reproduce the existing data by superposition of 
orbits. This means that restricting the viewing angle from 
the stellar kinematics alone will be challenging. 

To meet this challenge, more kinematic data in the in- 
ner Galaxy are urgently needed. It may well be possible to 
extract more information by studying the variation of kine- 
matic quantities with distance along the line of sight using 
well-defined selection functions. For external galaxies, where 
all stars are at roughly the same distance from the observer, 
it makes sense to record a single number for a kinematic 
quantity like the velocity dispersion along the line of sight. 
In studies of the Milky Way, this approach fails to exploit 
the full richness of information that is available to us. It 
would be more fruitful to calculate the kinematic quanti- 
ties for stars in the sample with, for example, magnitudes 
between m max and m m i n = m max — 1. As we have shown for 
the dynamical model, there are interesting and useful vari- 
ations of kinematic quantities with m max . These will prove 
invaluable in elucidating the structure of the bar, since the 
photometry together with the available line-of-sight disper- 
sions and streaming velocities are not enough to prescribe 
the bar uniquely. 
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APPENDIX A: NUMERICAL ALGORITHM 
FOR THE POTENTIAL 

This appendix sketches the calculation of the potential for 
the Binney et al. (1997) bar density. A straightforward ex- 



pansion in spherical harmonics converges rather slowly, es- 
pecially in the flat disk component. However, this problem 
can be overcome by a technique due to Kuijken & Dubinski 
(1994). Given a disk with density p — f(R)h(z), the potential 
is written as $ = $ + f(r)H(z), where r denotes spherical 
radius and H(z) is uniquely determined by H"(z) = h(z) 
and H(0) — H'(0) = (primes denoting derivatives). Insert- 
ing this into Poisson's equation gives 



V 2 $ 
4ttG 



[f(R)-f(r)]h(z) -f"(r)H(z) 



2 JM [H(z) + zH i {z)] . 



(Al) 



The expression on the right hand side is zero at z~0. The 
mass density generating <3> is not strongly flattened, so that 
it can be economically evaluated by, e.g., multipole expan- 
sion. Here, this technique is used for the two sub-disks. The 
subsequently evaluated potential $ and the modified den- 
sity is then expanded in spherical harmonics up to order 
I max = 64. The sharp truncation of this expansion at i max 
causes unphysical ringing of the resulting density and in or- 
der suppress this we have tapered the density expansion by 
multiplying it with exp( — [Z/32] 2 ). Finally, the potential and 
resulting forces have been computed on a pseudo-Cartesian 
101x101x161 grid of size 20 kpc x 20 kpc x 18 kpc. The points 
are linear in 2 ln(l+x/2), 2 In(l+j//2), and ln(l + 2z)/2. Po- 
tential and forces are then evaluated via a fifth order three 
dimensional spline. At each grid point, this spline gives the 
stored values of potential and its derivatives. The forces have 
everywhere continuous first and second derivatives. In par- 
ticular, the interpolated forces agree identically with deriva- 
tives of the interpolated potential, as is necessary if the Ja- 
cobi energy is to be accurately conserved along numerically 
integrated orbits. The evaluation of potential and forces is 
quick once the spline coefficients have been pre-computed. 



APPENDIX B: CALCULATION OF THE 
SECOND VELOCITY MOMENTS 

This appendix briefly sketches how to calculate the second 
moments for the regular and isotropic parts of the DF in 
turn. 



Bl Regular part of DF 

To compute a velocity moment projected along the lin e of 



2.4, we 



(Bl) 



sight towards (£o, bo) with the formalism of Section 
define 

n(to) = 8(£ - to) 5(b - bo) e(«) g(w). 

where g(w) is determined by the velocity moment in ques- 
tion. For example, to obtain «i os , one takes g — p-s, where s 
is the unit vector in direction of the line of sight. With this 
choice, one obtains for the projected moment using equa- 
tions £h to (|Tll) 



n[T 



1 



dte( Si (t)) g(wi(t)) 



x5(£i(t)-£ )5(bi(t)-b ) 



Here, 5(x) denotes a "<5-function" of some finite width, which 
in Baade's Window, for example, is 30'. The reason that we 



© 0000 RAS, MNRAS 000, 000-000 



A Dynamical Model of the Inner Galaxy 19 



cannot use strict 5-functions is that the probability for any 
single orbit to hit the line of sight in our finite integration 
interval of 200 dynamical times is zero. That is, the need for 
the finite-width functions 8 arises because we replaced (i) 
the integral over action-space by a finite sum (by virtue of 
the Monte-Carlo method), and (ii) the integral over angles 
by a finite, rather than infinite, time integral. 



B2 Isotropic part of DF 

Computing moments of the isotropic part of the DF / ls ° can 
be done directly since / lso is explicitly given as a function 
of Jacobi energy Ej and therefore phase-space coordinates 
w. The projected luminosity density is 



v(l,b)=Aix \ dss 2 €(s) / dvv 2 f(r(s,£,b),v) 



(B3) 



The second integral can be performed analytically for the 
second order splines Bi(Ej) that constitute the building 
blocks of / lso - see equation (^). The integration over s 
is done numerically. 

Since Ej depends only on the modulus of the velocity 
(and not v itself), the first moments J dvvf' BO vanish and 
one obtains (vi) = — v@i, i.e. the solar motions, while 



(B4) 



which is just the solar reflex motion VQi/s weighted by the 
density along the line of sight. For the second velocity mo- 
ments one finds (p= |p|) 



d SS 2 e (s) / dpp 4 r°(r(s,£,b),p), (B5) 



■In 



where the inner integral can be performed analytically by 
the aid of 



dpp A Bi(Ej) = 



35 AEj 



[Dti ~ 2DJ + Dj +1 ] 



(B6) 



with D% = Eji — $ c ff(r). The outer integral is performed 
numerically. For dispersions involving proper motions rather 
than velocities, one has to change the integrand by factors 
of a -1 . 



i.e. g is the density of states at given Ej. With 
d 3 p = 2irR~ 1 dE dL z (after integrating out one angle in p- 
space) and Ej = H — flL z (equation [15|) (Cf ) becomes 



g(E a \Ej) 



2tt 
IT 



2tt 

IT 



d 3 r 



4>(oo) 



dE 8{E-E a ) 



#(r) 



dL z 8 ( L z — 



Ry/2{E-<£) 



d0 



(C3) 



dz (i? 2 - />'i ) 

Here, Ri/ 2 (Eo,Ej,z,cj>) are the roots of 
(E - E 3 f = 2R 2 Q 2 (E - z, (/>)), 

more precisely, the minimum of these roots and the corota- 
tion radius - recall that the isotropic component lives only 
inside corotation. 



(C4) 



APPENDIX C: DIFFERENTIAL MASS IN THE 
ISOTROPIC COMPONENT 

We want to calculate the differential mass dM/dE, where E 
denotes the energy, for the isotropic part f lso (Ej) of the dis- 
tribution function. Let us consider a J-function at _Ej = Ej , 
then its contribution to dM/dE at E — Eo is simply the vol- 
ume of the four-dimensional cross-section of the phase-space 
hyper-planes at Ej = Ej and E — Eq 

g[E \E 30 ) = J d c 'w 5(Ej - Ej ) 8{E - E ). (Cf ) 



Then 

dM 
dE 



^ = Ma I dE, g(E\Ej) r°(Ej), (C2) 
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